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Abstract 

It is shown how to obtain accurate values for American options us- 
ing Monte Carlo simulation. The main feature of the novel algorithm 
consists of tracking the boundary between exercise and hold regions 
via optimization of a certain payoff function. We compare estimates 
from simulation for some types of claims with results from binomial 
tree calculations and find very good agreement. The novel method 
allows to calculate so far untractable path-dependent option values. 
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Pricing of options is an important area of research in the finance commu- 
nity. The field has been pioneered by Black and Scholes (1973). They made 
a major breakthrough by deriving a formula for the price of any contingent 
claim which depends on a non-dividend-paying stock. Using the assump- 
tion of no arbitrage, they were able to show that the price of a derivative 
security can be expressed as the expected value of its discounted payoffs. 
The expectation is taken under the assumption of a risk-neutral evolution of 
the value of the underlying security. Merton (1973) generalized these ideas 
to situations in which the interest rates themselves fluctuate in time. Due 
to the complexity of the underlying dynamics, numerical methods have be- 
come increasingly popular in modern finance. They are used for a variety 
of purposes, for instance valuation of securities and stress testing of portfo- 
lios. Analytical solutions for problems in finance have been found only for 
rather special cases. For instance, in order to gain the solution for prices 
of European options Black and Scholes had to assume that the evolution of 
the underlying asset price (more precisely its natural logarithm) follows a 
so-called Wiener process with time-independent volatility. European options 
can be exercised only at expiration date. In contrast, American-style options 
which can be exercised during the whole lifetime of the option can be valued 
only numerically. Furthermore, numerical methods have to be employed if 
the dimension of the problem increases (more than one state variable) or 
in case of more realistic approaches to the stochastic process like for some 
no-arbitrage models of interest rate evolution. 

The Monte Carlo approach lends itself very naturally to the evaluation of 
security prices and interest rates. Schematically, it consists of the following 
steps. First, sample paths of the state variables (asset prices and so on) over 
the relevant time horizon are generated. The cash flows of the securities on 
each sample path are evaluated, as determined by the structure of the se- 
curity in question. The discounted cash flows are averaged over the sample 
paths. The Monte Carlo method is very flexible, since it does not depend 
much on the specific nature of the underlying stochastic process. Its accu- 
racy is also independent on the dimensionality of the problem which is its 
dominant advantage over more traditional numerical integration methods. 
It is outside the scope of this paper to discuss the various facets of recent 
research on the use of Monte Carlo in the finance area. The reader may con- 
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suit the concise recent reviews of Boyle, Broadie, and Glassermann (1997) on 
these topics. In particular, they report on recent progress in developing more 
efficient Monte Carlo algorithms. Standard Monte Carlo methods converge 
notoriously slow (with 1/y/N, N being the number of sampled paths). A 
central goal of recent research activity has therefore been the refinement of 
variance reduction techniques (antithetic and control variates), importance 
sampling and low-discrepancy random number sequences. If the development 
in sciences like physics is any guide, Monte Carlo applications in finance will 
become even more important. With computational prowess further increas- 
ing the Monte Carlo method may even shed its 'brute force' image from the 
past which was founded on its rather slow convergence property. 

Boyle (1977) was the first to propose the use of Monte Carlo methods 
for option pricing in the literature. Since the exercise date of European 
contingent claims is fixed, the mechanics of a price evaluation of European 
options employing the Monte Carlo method is rather straightforward. Let us 
begin by collecting the notation for the most important variables which we 
are going to use throughout the rest of the paper. In general, we follow the 
notation in Hull's book (1993), except that the present time for which the 
option price is calculated is designated to be "0" (without loss of generality). 

S t = the value of the state (asset price, . . . ) at time t 

underlying the derivative security 
(which may be a single variable or a vector), 

So = the initial value of the state variables, 

X = strike price of the option, 

p = path on which the underlying state evolves in time, 

Gtip) — the value of the option if exercised at time t 
(max(S t — X, 0) for a simple Call and 
max(X - S t ,0) for a Put), 

r t = instantaneous risk-less interest rate at time t, 

a = volatility (standard deviation) of state variable 

or their components respectively, 
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T = the lifetime of the option contract (expiration date), 

Vt = the option value at time t. 

We may easily estimate the European option value from a Monte Carlo 
path sample for the variables related to the derivative security. These paths 
need to be generated with probabilities determined by the underlying stochas- 
tic process. The option value at the expiration date T has been specified in 
the option contract. In general, it may depend on the path of the relevant 
state variables between closing the contract and expiration date. Each path- 
wise final option value from the simulation is discounted backward in time 
to determine its value at time 0. Since we do not know at time t=0 yet on 
which path the underlying variables will evolve, the average over all paths is 
taken to estimate the present option value: 

N 

K,^iV" 1 ^exp(-rT)\/ T (p) (1) 
P =i 

The discount factor includes the risk-free interest rate r, possibly averaged 
over the lifetime T of the option. For simplicity, it is assumed in eq. (|1|) 
that all sampled paths have equal probability. Generalization to a situation 
in which paths are characterized by different probabilities which may arise 
e.g. in the context of Monte Carlo optimization by importance sampling is 
straightforward. 

The pricing of options which may be exercised prior to the expiration 
date by Monte Carlo simulation is more involved. Here, the owner holds the 
right to exercise the option at several (Bermudan) or possibly infinitely many 
(American) 'decision dates'. Many types of American contingent claims trade 
on exchanges and in the over-the-counter market. Examples include options, 
swaptions, binary options and Asian options. It has been shown by Roll 
(1977), Geske (1979), and Whaley (1981) that exercise of call options which 
give the owner the right to purchase some underlying asset at the agreed-on 
strike price is usually unfavorable before expiration date, except close to ex- 
dividend dates. The situation is much less clear-cut in case of put options 
which give the owner the right to sell the asset at the strike price. In these 
situations the 'naive' Monte Carlo is encountering unsurmountable difficul- 
ties. In order to decide whether it is more favorable at some intermediate 
time to hold or exercise the option the owner needs to compare the expected 
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payoff in the two cases. The maximum of the two forms the option value at 
that time which - after discounting and taking the expectation value - leads 
to the following expression for the option at closing date: 

V = (exp(-rT)VV) . (2) 

The angular brackets denote the expectation value with respect to the (risk- 
free) probability measure. The path-dependent times T are the so-called 
optimum stopping times which may be any of the decision dates. 

A single path provides clearly insufficient information to evaluate the 
option value in case of non-exercise at any of the decision dates. The insuffi- 
ciency of the naive Monte Carlo method to deal with the optimum stopping 
problem has lead some authors like Hull (1993) to the claim in the literature 
that "Monte Carlo simulation can only be used for European-style options" . 
On the other side, evaluation of American-style Options via Monte Carlo sim- 
ulation has found already some consideration in the literature as reviewed 
Boyle, Broadie, and Glassermann (1997). Their common denominator is 
a 'clever' estimate of the option prices on the decision dates by bundling 
subsets of paths. These authors noted that the suggested algorithms can- 
not be considered satisfactory yet. Either some approximations may lead to 
uncontroled errors affecting the simulation results or the required large com- 
putational effort limits their applicability, e.g. to a small number of exercise 
dates. 

In this paper we are suggesting a completely different strategy to calcu- 
late American-style put and call option prices via Monte Carlo simulation. 
During the sampling of the paths we do not attempt to estimate the expec- 
tation of the payoff if the owner continues to hold the option. Instead, we 
use the sampled paths to evaluate the boundary between the early-exercise 
region and the hold region in the space of variables entering into the option 
contract. The location of this boundary is the crucial piece of information 
whose knowledge allows the straightforward use of the Monte Carlo proce- 
dure for option price estimation. We treat the valuation of American options 
as an optimization problem of a certain payoff function which depends on 
the set of sampled paths. This function depends also on the exercise policy, 
i.e. the boundary between the early-exercise region and the complementary 
hold region (YES-NO boundary in the following). Maximation of the payoff 
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function provides an estimator of the YES-NO boundary which gets arbi- 
trarily close to the true location of the boundary with increasing number of 
sampled paths. For the plain vanilla put and call options, the boundary is 
rather simple at any exercise date, a point for one state variable, a line for 
a two-dimensional space of state variables and so on. After the boundary 
has been estimated the price estimation proceeds as for European options, 
because the option prices for points on the boundary are known. The only 
difference to price estimation of European options as in eq. ([!]) is that - con- 
cerning paths crossing the YES-NO boundary - each pathwise option value 
is discounted starting from the point at which the path crosses the boundary 
for the first time: 

N 

V^N- l ^M-rr{p))GrUp) ■ (3) 
P =i 

Time T(p) denotes the earliest time at which the p-th path crosses the YES- 
NO boundary (if at all). Otherwise T(p) equals the time until expiration of 
the option (T) . Gr(p) (p) is the agreed-on payoff from exercising the option at 
time T(p) under the assumption that path p represents the realized evolution 
of the underlying state variables. Here we have tacitly assumed that the 
starting point lies in the NO region. Of course, the algorithm is able to 
handle the other possibility - immediate exercise - as well. In this case 
calculation of the option price is trivial, however. 

In a later section of this paper we will demonstrate that the novel Monte 
Carlo algorithm achieves an accuracy in the determination of American op- 
tion values which compares well with corresponding calculations for Euro- 
pean options. We compare the results for a spectrum of standard American 
Puts with results from binomial tree calculations. This method originally 
introduced by Cox, Ross, Rubinstein (1979) is widely used for option valu- 
ation. Here we assume an ideal Wiener-type stochastic process for a single 
state variable, because a binomial tree is well suited to provide very accurate 
results in this situation. Of course, this comparison serves mostly illustra- 
tive purposes to show that the method works at finite N. Lateron, we are 
going to present numerical results for path- dependent American options. We 
consider two cases, Puts on either the geometric or the arithmetic average 
over the lifetime of the option. Assuming again an underlying ideal Wiener 
process we may compare the Monte Carlo estimations to results from bino- 
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mial tree calculations. While options on the geometric average can be priced 
by the tree method to arbitrary accuracy no such tree method is available 
for options contingent on the arithmetic average value. Some approximative 
method has been suggested to express the latter option prices by the corre- 
sponding prices in case of geometric averaging. Therefore we will compare 
the approximative formula to our more accurate procedure. In general, it is 
envisioned that the novel Monte Carlo approach and some variants will be 
applied for valuation problems which cannot or not easily be tackled by other 
methods. Such problems encompass stochastic processes with Non-Markov 
properties, stochastic volatilities, multi-dimensional state variables and other 
path-dependent American options. 

I Valuation of American- style options as an 
optimization problem 

Let us assume in the following that the underlying state variables whose 
values determine the value of an option at exercise dates follow a known 
stochastic process. Generally, a stochastic process is specified by defining 
the state space, an index parameter (usually the time) and the dependence 
relation between the random variables. The latter is usually given in terms 
of a stochastic integro-differential equation containing some deterministic 
drift terms and random changes in addition, notably diffusion or jump pro- 
cesses. Memory or so-called non-Markov effects may readily be included into 
the transition probabilities. Non-Markov models have found recently con- 
siderable interest in the finance community, because the Heath, Jarrow and 
Morton (1992) approach to interest rate evolution can be shown to possess 
Non-Markovian features, in general. We assume 'smoothness' for the transi- 
tion probability density function which determines the probabilities - based 
on a path's evolution history - that the path will comes close to any element 
of the state space at some time in the future. These are supposed to be con- 
tinuous functions of the path variables. A discrete set of sequentially selected 
random points in time is often called a "chain" . A path is the corresponding 
generalization for a continuous time evolution. 

At first, we will restrict ourselves to options and stochastic processes 
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which result in the following simple property of the YES-NO boundary: 

Al: "The boundary between YES and NO regions at any given decision 
time t < T is a simply connected hypersurface in the state space with 
dimension of one unit less than the dimension of the state space." 

As a consequence of the assumption, there are no disconnected pieces 
of YES or NO regions at a given time either. With the exception of some 
pathological cases, assumption Al will be compatible only with option payoffs 
contingent on the values of the state variables at exercise date. Instead of 
specifying these pathological stochastic processes we simply require here: 

A2: "The payoff from exercising an option is only a function of the state 
variables at the exercise date but not at earlier times, i.e. G t {p) = 
G t (S t ,X)." 

In general, the topology of exercise and hold regions of American options 
may be more complicated than allowed by our assumptions. Two examples 
are some types of barrier options and path-dependent options. Suppose the 
owner of a barrier option receives only some payoff if a security stays within 
some specified range of values at exercise date. Therefore the option will not 
be exercised for 'too small' and 'too large' values of the underlying which leads 
to two disconnected NO regions. Another perhaps more interesting example 
of options not covered by the assumption stated above are American path- 
dependent options whose valuation we will return to later in this section. 
Here the payoff depends on the values which the underlying state variables 
have taken in the past, e.g. on the arithmetic mean within some time interval. 
Projected into the space of state variables at any given time YES and NO 
region are completely scrambled. It may be optimal to exercise or hold an 
option for the same present value of the state variables S t depending on the 
different histories of each path. Lateron, we will show that a redefinition of 
the state space enables us to cover the case of path-dependent options by the 
same algorithm. 

We should also mention that the assumption A2 restricts the memory 
properties of the underlying stochastic process. They should not spoil the 
assumed property of the YES-NO boundary in the state space. This simply 
means that although the future evolution of each path may depend on the 
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past, the 'average outcome' determining the decision whether to exercise the 
option at present time is not influenced by the past but of Markovian nature. 

Assumption Al allows us to define coordinates in the ds dimensional 
state space such that any element St of the state space can be uniquely 
characterized by some point of the YES-NO boundary Sf and an additional 
coordinate c t 

St = (S?,c t ) (4) 

in a peculiar way. For any path passing through St at present time the 
following inequalities hold between the payoff from exercising the option G t 
and the option value V t N if the option is being held at least until the next 
exercise date: 

( > V t N for ct < 
G t {S t ,X)\ <V t N for Q >0 (5) 
[ = V t N for c t = 

The content of eq. @ is to provide convenient coordinates from which one 
can easily read off whether point St is a boundary point (c t = 0), in the YES 
(q < 0) or in the NO (c t > 0) region. The knowledge of all c t or just their 
signs is tantamount to finding the YES-NO boundary itself. Unfortunately, 
a straightforward application of the inequalities as expressed in eq. (0) to 
determine the boundary is prohibitive. Estimating the expected payoff from 
holding the option and comparing it to the payoff from exercising would re- 
quire a 'new simulation within the simulation'. Such a procedure becomes 
numerically awkward and would severely restrict the complexity of the prob- 
lem to be tackled, concerning e.g. the nature of the stochastic process and 
the number of exercise dates. 

We will proceed now to estimate the location of the YES-NO boundary 
based on a set of sampled paths. For this purpose it will be necessary to 
evolve the state variables backward in time. This feature of the novel Monte 
Carlo approach which it shares with all other numerical approaches to option 
pricing like the binomial tree method is caused by the boundary conditions 
of the problem. The option price is a specified function of the state variables 
only on the expiration date. Of course, the location of the YES-NO boundary 
is also trivially determined at expiration time. For a simple call or put option 
on a single asset it is the point St = X. 

The continuous time process will be presented as a sequence of time steps 
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Sti, i — 1, . . . Nt which add up to the lifetime of the contract J2 St% = T. The 
time steps are supposed to match the distance between decision times exactly 
if a finite number of exercise opportunities has been specified in the option 
contract. For American options which give the owner the right of exercise at 
any time this is an approximation. However, choosing 'small' time steps St 
the results will be close to the continuum limit. The current time t for which 
the location of the YES-NO boundary is being sought can be represented by 
the integer index i (with t = i-5t for American options). Since the location of 
the YES-NO boundary will be determined step by step moving backward in 
time, it is implicitly assumed from now on that the YES-NO boundary has 
already been determined for later times z+l,z+2,. . . ,Nt (within the usual 
uncertainties due to the numerical imprecision). Therefore we are able to 
evaluate for each path 2 its payoff from optimal exercise at later times: 



Gi+iip 



>t) 



G t (S i+ i(p),X) f° r S i+1 (p) in YES region 

exp(-££=i + ir fc <5t fe )- 

GtiS&lX) for Si+i(p) in NO region . 

(6) 

The path-dependent payoff has been discounted to time i + 1 in eq. (||. j 
represents the earliest time (under the constraint j > i) at which path p 
crosses the boundary from the NO region into the YES region or equals Nt 
if the path stays within the NO region for all times between t and T. (For 
the latter case G t (Sj(p),X)=0, of course.) We have used the subscript > t 
for the argument p of Qi to emphasize that the payoff defined above depends 
only on the values of the path variables for times > t but not on the past 
evolution of the state. So far, we are not able to calculate the discounted 
pathwise payoff Gi(p>t), i-e. for time t . This would require the knowledge 
of the YES-NO boundary at time t. However, the path-dependent payoffs 
given in eq. (|6]) enter into the expression for the value of the option at time 
i subject to the condition that the option is not immediately exercised (and, 
of course, that the state variables have evolved along path p so far): 



Vf = exp(-r^) J Dp >t (7) 



2 For the ease of reading we are using calligraphic letters in this paper when dealing 
with functions of paths or path segments. 
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Gi+i(p>t) 7r(p >t ;p< t ,t) 



The integration over all future path segments in eq. (0) 




is an ordinary multiple integral (but a path integral in the continuous time 
limit). 7i(p >t ;p< t ,t) denotes the transition probability density function and 
describes the probability (density) to evolve along the path segment p >t in 
the future up to time T depending on the present - and previous (for Non- 
Markov processes) - states of path p< t . 

Next we define a continuous path H in the state space by letting c be 
a continuous parameter within some range such that all S'(c) = (Sj , c) are 
elements of the state space and vice versa. Let us assume for a moment 
that we know already that this line intersects with the YES-NO boundary 
but we do not know which point it is. (Somebody else has constructed the 
path H.) Can we estimate the value of the YES-NO boundary point using 
the information provided by the path sample and obtain its correct value 
in the limit iV — > oo? The answer is positive. We are able to construct a 
function of the states S'(c) e H whose only maximum is at S'(c) = (S]f,0), 
i.e. the boundary point, for iV — ► oo. It is clear that our ability to define 
such a function has important repercussions for the question how to find the 
YES-NO boundary. The virtue of path H is to provide us with the knowledge 
that it has exactly one intersection point with the YES-NO boundary. We do 
not necessarily need this information. Instead, we may search for the global 
maximum on continuous paths connecting points in the YES region with 
points in the NO region. According to statement Al they will cut through 
the YES-NO boundary at least once. (It will not hurt the search procedure, 
if we accidentally hit the boundary several times reflected by several local 
maxima.) It is also rather straightforward to determine selected points which 
are in either the YES or the NO region for sure. For instance, points lie in 
the NO region if the payoff from exercising is zero. Or the point with zero 
value of the underlying is certainly an element of the YES region for a plain 
vanilla Put, because the payoff cannot get any better. 

In order to proceed with the proof of our conjecture we define a - generally 
suboptimal - exercise policy for all paths intersecting with H at time t which 
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depends on an arbitrary point S' = (5* t y , c') G H. We define a preliminary 
"YES" region (with the quotation marks reminding us that the true YES 
boundary may be different) as to cover all paths with state variable Si(p) = 
(SY,Ci(p)) which fulfill Clip) < d. For instance, in case of a plain vanilla 
put (call) option on a single underlying asset we may declare that Si(p) lies 
in the 'YES" region if it fulfills the condition Si(p) <(>) S'. The pathwise 
payoff from the 5"-dependent policy is 

rry ( ai\ \ G t(Si{p),X) for Ci (p)<d, 

P - fe " S) - = lexp(- rA )ft +1 (p>«) else . (8) 

It is clear that the pathwise payoffs will sometimes be larger for 'wrong' 
boundary points {d 7^ 0) than for the correct one. This reflects the additional 
information about the future which is encoded in the path segment p> t - The 
essential point is now to sum over all paths which removes this bias: 

N(H) 

0?(S'):=N(H)-' £ V t (p> t ,S>) . (9) 

v 

The function Of(S') depends on time, the sampled paths intersecting with 
H (of number N(H)) and on S'. Our central proposition can be stated now 
as follows: 

PI: "Taking the limit of infinitely many paths, 0^°(S') := lim^^H)^oo 

has one extremum which is a maximum within the one-dimensional 
domain H . Its position is the intersection point of H with the YES- 
NO boundary S'=(Sf,0)." 

In order to prove proposition PI we are going to rearrange the sum in 
eq. @. We put all paths which have essentially the same transition proba- 
bility density function into the same (albeit infinitesimally small) 'bin'. Here 
we need the assumption stated above that the transition probability density 
function is a continuous function of the path variables. After 'binning' and 
taking the limit N(H) 00 eq. (P) turns into 

Or(S') = J Dp< t p(p< t ) (J Dp >t V r (p> t ,S')n(p >t ;p< t ,t^ . (10) 



The transition probability density function ir(p >t ;p<t,t) appears in eq. (|T0|), 
because relative frequencies approach the probabilities which are governing 
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the random selection of paths according to Bernoulli's law of large numbers. 
The probability density p(p<t) to have a path segment close to p< t needs to be 
defined in accordance with the integration measure Dp< t and is normalized 
to one. 

The quantity in brackets () on the right hand side of eq. ( |TD| ) 

D t (S,S') := J Dp^V^S^^p^p^t) (11) 
can be readily evaluated using eqs. (|7D and (§) to be 



D t (S,S') = { °i S > X) fOTC<C ' (12) 
^ V t N {S) else . 



D t (S, 5") coincides with the option price V t (S), except for c values in a certain 
range. Furthermore, on virtue of eq. (|5]) D t (S,S') is always smaller than 
V t (S) for such chosen c values: 



G t {S, X) for c G (0, d) with d > 0, 
Vf(S) c G (d, 0) with d < 
D t (S, S') else 



The conditions of eq. (^) together with the positiveness of the weight factors 



p in eq. (|ToD guarantee that 

0?(c a ) < 0™(c b ) hrc a <c b <0 (14) 
O t °°(c ) < O t °°(c b ) forc a >c 6 >0 . 

This completes the proof of proposition PI. 

As stated above we may turn things around and construct the YES-NO 
boundary from a set of lines in state space which connect elements of the 
YES and NO region. The global maximum of function O^(S') along each 
line is the supposed intersection point with the YES-NO boundary. This 
process may be iterated until the desired precision is reached. This standard 
optimization procedure becomes particularly simple in case of plain vanilla 
options. Since the YES-NO boundary is just a point at any exercise date, 
its location is already completely specified with one iteration step. The 
determination of the YES-NO boundary at time t (index i) enables us to 
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calculate the pathwise payoff at this time. Going step-by-step backward 
in time we may finally calculate the pathwise payoffs at initial time whose 
average determines the option value at contract time. 

The presented optimization technique may be used in the problem of val- 
uation of American path- dependent options as well. Here the payoff from ex- 
ercising the option depends on some function of the state variables in the past 
which we will denote as St(p<) in the following. For instance, St(p<) may be 
the - arithmetic or geometric - average of any state variable over some time 
window or the value at some specified time(s) of the past (look-back options). 
For simplicity, we restrict ourselves here to stochastic processes of Markovian 
nature whose parameters have been specified. In this case the exercise de- 
cision at any time before expiry is a function of just two - each possibly 
vector-valued - variables, the values of the state variables and of the path- 
dependent variable St(p<), both taken at the decision time. The difference 
between standard path-dependent and plain vanilla options is just that the 
role of the current value of the underlying variable in plain- vanilla contracts 
is played by S t (p<) in path-dependent contracts, e.g. Max(X — S t (p<), 0) for 
a Put. The topology of the exercise region for this kind of path-dependent 
options is very simple but only in the space of state variables enlarged by the 
additional dimension(s) of S t {p<) ■ The YES-NO boundary in the S t -St{p<) 
space is again simply connected. Only its projection into the St space leads 
to a scrambling of YES and NO regions. This suggests immediately how 
the suggested optimization technique may be applied for the case of these 
path-dependent options. We equip the space of states with additional dimen- 
sions by declaring St(p<) to be one of the stochastic variables. Of course, 
the stochastic character of this variable is rather funny. Its change with time 
is completely deterministic and of Non-Markovian type. However, these fea- 
tures are in accordance with the basic assumptions entering into the proof of 
proposition PI. The complications of such path-dependent option contracts 
just result in a higher dimension of the space in which we have to track the 
YES-NO boundary. 

We are now in a position to compare our novel strategy for an extraction of 
American-style option prices from Monte Carlo simulation to previously sug- 
gested ones. Boyle, Broadie, and Glassermann discuss mainly three strategies 
in their review paper to overcome the inherent limitations of the naive Monte 
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Carlo procedure: the bundling algorithm invented by Tilley (1993), the strat- 
ified state aggregation algorithm suggested by Barraquand and Martineau 
(1995) and Broadie and Glassermann's (1995) algorithm based on simulated 
trees. In Tilley 's approach bundles of paths which are 'close' in state space at 
given time are considered to emanate from one single point for the purpose 
of option valuation. This introduces some error whose magnitude cannot 
be easily controled. The bias of the exercise decision will be very strong if 
a bundle contains only few paths. Such a situation is likely to emerge in 
situations in which there are many exercise dates or several state variables. 
Barraquand and Martineau introduce another type of 'bundling', however, 
in the payoff space. Our approach shares similarities with Tilley's algorithm 
in the sense that it contains an implicit 'bundling' of paths (called 'binning') 
in the proof of the proposition PI. However, we have actually never to carry 
out this bundling in a numerical calculation, because we stick to an opti- 
mization of the 'inclusive' function 0%°(S') which sums up all bins. Broadie 
and Glassermann proposed an algorithm based on a 'bushy tree' structure 
in state space which avoids partitioning of the state or payoff space. In their 
model many branches emanate from each node. This process is replicated as 
often as there are decision times d in the problem. The replication process 
effectively limits d to small numbers (on the order of 4) in practical compu- 
tations. They also discuss the problem that the suggested estimators of the 
option prices tend to be biased in 'up' or 'down' direction, and converge only 
to the correct value in the limit N — > oo. This problem is of great relevance 
for practical calculations (which are always at finite N). They use the idea 
to calculate two estimators - one biased low and the other high - to obtain 
confidence limits from the simulations. Lateron, we are going to apply this 
idea in the framework of our approach. 



II Simulation results for American options 

In this paper we will only consider stochastic processes of rather simple na- 
ture, so-called ideal Wiener or diffusion processes (for the log of any state 
variables). In physics these processes are usually called Brownian motion. 
Such processes seem to be most commonly used in valuation problems. Fur- 
thermore, numerical techniques have been developed to calculate prices for 
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plain vanilla and other simple types of American options. In particular, the 
binomial tree method may provide very accurate values and will be used as a 
benchmark here against which the results from the Monte Carlo simulations 
are being tested. In this section we will restrict ourselves to Puts on the 
current value and on the - geometric or arithmetic - average of one single 
underlying security. An ideal Wiener process can be characterized by the 
differential transition law 



dt > 0. z is a random number drawn from a normal distribution with mean 
and variance 1. /i is in general the drift velocity which may be replaced by the 
risk-free instantaneous interest rate r using arbitrage arguments. Sometimes 
one needs to follow an ideal Wiener process only backward in time, e.g. in 
case of Monte Carlo estimation of plain vanilla options. Using the bridge 
construction - start So and final value St are known - a random value at 
intermediate time may be chosen according to 



After having fixed final and start values, iterative use of eq. (|16D may be 
employed to generate a time reverted copy of the standard diffusive Wiener- 
type motion. For the simulation results presented below we have used a 
simple importance sampling scheme. The exercise value of strongly out- 
of-the money options at expiry will be zero on most paths. Therefore the 
weight of the rare paths resulting in nonzero option value has been artificially 
enhanced in the path generation process. Of course, the enhancement factors 
are taken out again in the calculation of the path-averaged prices. 

In its most simple form the Monte Carlo algorithm which is suitable for 
valuation of American-style options consists of the following steps: 

(1) Generate path sample. 

(2) Go step-by-step backward in time starting at expiry. 

(3) Track the YES-NO boundary by employing the optimization of func- 

tion Of(S') along arbitrary paths connecting YES and NO region (cf. 




(15) 



governing the change of the state variable St := In St within a time interval 




(16) 



eq. ©). 



16 



(4) Evaluate the pathwise payoffs at the earlier time analogous to eq. (||). 

(5) Repeat the process until the contract time is reached. 

(6) Average the pathwise payoffs discounted to initial time over the whole 

path sample (cf. eq. (||)). 

Let us discuss now some aspects related to the finite size of the path sam- 
ple in practical calculations. It is clear that step (6) does not provide an 
unbiased estimator of the option price at any finite N but is biased upward. 
The reason is that an error is made extracting the YES-NO boundary from 
the finite size path sample. A deviation of the extracted from the true bound- 
ary has its root in the information about the mismatch between frequency 
distribution of the N paths and underlying probabilities which enters into 
the optimization process. The estimator becomes only asymptotically unbi- 
ased. Since we are not able to construct an unbiased estimator, we may get 
an estimate of the error introduced by the bias in step (6) by constructing 
an estimator which is downward biased at finite N but also asymptotically 
unbiased. It is very simple to construct such an estimator: 

(6') Generate a path sample independent from the first one and evaluate 
the American option prices using the previously calculated YES-NO 
boundary. 

The reason that the prices calculated from this path sample are biased down- 
ward mirrors the one for the upward bias in the former case. Here any 
deviation from the true boundary implies a suboptimal exercise policy. 

We hasten to add that the algorithm presented so far has one shortcoming 
related to the combined effects of finite sample size and well-determined 
values of the underlying variables at start time. Only with paths crossing 
the YES-NO boundary its optimal location can be determined, of course. 
Formally, the proof of proposition PI requires that the probability densities 
appearing in eq. (|10|) be positive-definite in a region covering the YES-NO 
boundary. However, this condition will be violated in practice at early times, 
since all sampled paths emanate from a common starting point. On the 
other side, in such a situation the knowledge of the precise location of the 
boundary becomes irrelevant for the valuation problem. It suffices to know 
at this stage on which side of the boundary the vast majority of paths can 
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be found. Two variants with differing degree of sophistication have been 
developed to determine the earliest time for which the algorithm above is to 
be used: 

(7a) If the YES-NO boundary gets closest to either one of the endpoints 
of H, e.g. the paths with smallest or largest St value for plain vanilla 
options, it is assumed that the boundary can be found outside the range 
of the path sample for all earlier times. 

(7b) "Flashlight mechanism": when the YES-NO boundary moves towards 
a region only poorly covered by the original path sample, additional 
path segments are created which evolve in the region of the YES-NO 
boundary at this time. As far as the proof of proposition PI is con- 
cerned the probability densities in eq. ( |10D do not need to be consistent 
with the initial value specification for the pricing problem. Here we 
make use of this freedom to generate a set of paths whose initial values 
are tuned to 'shed light' on the location of the YES-NO boundary. 

Concerning accuracy simulation mode (7b) is only barely noticably superior 
over mode (7a) for the tested sample of plain vanilla options. All simulation 
results for more complicated options reported on in this paper have been 
achieved employing mode (7b) . 

At finite N the values and therefore also the maximum of the payoff 
function O^(S') do not change between values of its argument which are 
taken by two neighboring paths. Furthermore, the precision with which the 
true maximum of function Of°(S') can be estimated from a finite size sample 
is a multiple (^> 1) of the typical distance between paths in the region along 
the YES-NO boundary. Fig. [I] shows the Monte Carlo estimation of the 
optimization function Of(S') for a plain vanilla put option on a security of 
(time-dependent) value S in order to illustrate the finite N effect. Selected 
O^(S') values are displayed at a fixed time but varying the size of the path 
sample. The true minimum is quite shallow for the selected time as can be 
seen from the figure. The shallowness is related to the strong time dependence 
of the location of the YES-NO boundary, because we consider a time shortly 
before expiry. The location of the maximum can be reliably extracted from 
the simulations only for rather large samples on the order of 10 5 paths. One 
may ask whether in view of the imprecisions associated with any finite sample 
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size it is worthwile to search for the exact maximum of Of(S'). In fact, a 
well-chosen 'smoothing' procedure which subtracts the white random noise 
may give more accurate results than taking the real maximum of Of(S') as a 
boundary point. For this paper we compared two variants how to determine 
the YES-NO boundary points: 

(3a) Determine the YES-NO boundary point as the element of the set of all 
intersection points Si(p) of the sampled paths with path H for which 
Of restricted to these points attains a maximum value max{s.( p )} Of . 

(3b) Determine the YES-NO boundary point as the element of the set of 
all grid points Sf for which Of restricted to these points attains a 
maximum value m&x^s^G)} Of . The distance between sites on the grid 
is lowered stepwise until some pre-defined precision is reached or more 
than one local extremum appears for the set of values of Of restricted 
to the grid sites. 

Ultimately, in the limit of infinitely many paths both variants may give the 
same correct result for the optimization of Of. However, computation em- 
ploying variant (3b) is much quicker. It requires a number of operations 
growing only linearly with N while variant (3a) needs at least on the order 
of iVlogiV operations. 

We have generated Monte Carlo estimates for prices of randomly sam- 
pled standard American-style put options (see Table I for details). We have 
employed the two simulation modes (3a) and (3b) to compare their ef- 
fectiveness. The frequency distribution dPj dx of the relative errors x in the 
simulation results based on Monte Carlo mode (3a) are shown in Fig. [| 
for mode (3b) in Fig. |3|. These errors are determined from the normalized 
difference to the prices extracted from binomial-tree calculations: 

x = (V» G -V*)/V* . 

The Monte Carlo estimates have been calculated for the path sample which 
has been used to determine the YES-NO boundary and for an independently 
generated path sample. The comparison between the results in the two modes 
reveals that both lead to similar accuracy. In addition, the figures display 
the distribution of the differences between Monte Carlo values based on the 
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original and the independent path sample. Upward and downward bias are 
clearly revealed for the simulation results in mode (3a) (lower left panel 
of Fig. 0). On the other side, the fact that mode (3b) employs only an 
approximate optimization of the payoff function Of tends to wash out most 
of the up- and downward biases. Fig. |3] shows also the error distribution 
of the corresponding European option estimations. We note that the typical 
size of the errors is comparable to the results for American puts. The error in 
both cases is thus dominated by the finite size of the path sample (N = 10 5 
in Figs. |2| and |3|). The error induced by the bias in the construction of the 
YES-NO boundary at finite N seems to be of less importance. 

Let us turn now to the case of put options contingent on the geometric 
or arithmetic mean value of the underlying variable. As discussed in the pre- 
ceeding section exercise and hold regions and therefore the YES-NO bound- 
ary are simply connected in the S t — S~t(p<) space. The boundary is therefore 
one-dimensional in this space. We have chosen a very simple procedure in 
order to track the boundary. These options are not directly contingent on 
the current S value at any of the decision times which makes St{p<) the more 
relevant variable. Therefore we sort the sampled paths into bins covering the 
relevant St region. Ignoring the differences between St values inside each bin 
we have reduced the optimization problem again to finding the St(p<) value 
at which the payoff function Of attains maximum value. This optimization 
proceeds as described already for plain vanilla options. For the calculations 
described here we have taken the number of bins in St direction to be 20. 

Monte Carlo simulation results for a randomly generated sample of put 
options contingent on geometric averages are presented in Fig. As before, 
errors in the simulation have been calculated in case of geometric averaging 
by comparing the Monte Carlo estimations for the option prices to binomial 
tree results. Here we consider only averages over the whole lifetime of the 
option. This is an n 3 problem for the binomial-tree method, n being the 
number of time steps. Allowing the time window for the averaging to slide 
would increase the storage and CPU time requirements by one power of n. 
In contrast, the Monte Carlo approach is not impacted severely by such a 
generalization. 

A comparison with the corresponding errors of the simulation results for 
European options under the same conditions (not shown) reveals again that 
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the errors are dominated by the finite size of the path sample and not by the 
uncertainty in determining the YES-NO boundary. The effect of upward and 
downward bias in the estimations due to choosing either the original or the 
independently generated path sample are very weak as can be seen from the 
lower left panel of Fig. [|. As explained before this is a direct result of the 
approximations in the optimization procedure (3b) . It is noteworthy that 
the bias from parametrizing the YES-NO boundary on a grid (in St direction) 
is stronger than either bias from taking one of the two path samples. The 
representation of the YES-NO boundary on a grid tends to underestimate 
the option prices because of the coarse graining procedure involved. The 
small nonstatistical downward net bias can be read off from the average 
deviation between Monte Carlo estimation and binomial tree results for the 
options contingent on geometric averages. It amounts to - 0.24 % for the 
results based on the original path sample and - 0.1 % for the independent 
path sample. 

Results for prices of options contingent on arithmetic averages over the 
whole lifetime of the option are also included in Fig. [| (lower right panel). 
Since no other accurate method is available, we present the distribution of 
relative deviations between the values based on an approximate formula sug- 
gested by Ritchken, Sankarasubramanian, and Vijh (1989) and the simulation 
results. The approximate formula relates the price for a claim contingent on 
the arithmetic average to the corresponding option price using geometric 
averaging via 

A am ~ Ji gm & gm -+- J^ am . { 1 ( ) 

Here A (E) denotes "value of American (European) option". The subscript 
am (gm) refers to arithmetic (geometric) averaging. The superscripts bt and 
MC characterize the method of calculation, binomial tree and Monte Carlo 
respectively. Fig. [| shows that the approximation works reasonably well. 

Ill Conclusions 

We have shown in this paper how to construct Monte Carlo algorithms for 
American option valuation. We are in the lucky situation that the suggested 
algorithm turns out to be conceptually simpler than the other algorithms 
suggested so far. It depends only linearly on the number of sampled paths 
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and exercise dates which makes it computationally feasable to get close to 
the continuum limit for these two variables. Thus no price is to be paid 
for the accuracy of this novel Monte Carlo algorithm in terms of more CPU 
time or exceedingly large storage requirements. The crucial element, the 
optimization of a certain payoff function of the path sample, turns out to be 
numerically rather stable. The tracking of the YES-NO boundary is achieved 
most precisely in those regions of the state space covered densely by sampled 
paths, i.e. where it is needed for the option price estimation. There is a 
further reason for the algorithm's stability. In case a maximum turns out 
to be shallow and is therefore not easily found by the algorithm, the option 
price does not depend much on the exact location of the boundary between 
early-exercise and hold region. Indeed, we have demonstrated that the errors 
in the simulation results are typically dominated by the statistical errors 
if the size of the path sample is on the order of 10 5 . This points to one 
important direction of future research. Paskov and Traub (1995) have shown 
that the replacement of pseudo-random numbers in the traditional Monte 
Carlo approach by a series of so-called quasi-random numbers may lead to a 
spectacular gain in the precision of the achievable results, at least for some 
problems. The reason is that the quasi-random numbers are distributed more 
evenly. By using appropriately chosen sequences of quasi-random numbers 
sample sizes for estimation of European options could be reduced by orders 
of magnitude, without loss of accuracy. It would open new dimensions for the 
applications of the novel Monte Carlo algorithm if this would also hold true 
for American options. Moreover, we would like to add that the extraction of 
sensitivity coefficients - usually called the "Greeks" - employing the novel 
Monte Carlo algorithm is as straightforward as for previously considered 
situations, e.g. by Fu and Hu (1995). 

We have presented some simulations for option types for which other 
means of calculation are known. However, the Monte Carlo algorithm pre- 
sented here provides unique opportunities to evaluate path-dependent option 
prices for which no other methods are available or computationally feasible. 
One example, the case of options contingent on arithmetic averages has been 
discussed already in this paper. Another candidate are "look-back" options 
of American style. Usually, look-back options allow the owner to buy or sell 
an asset for a price dependent on the values of the asset during the whole 
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lifetime until expiry, in the simplest case the minimum or maximum. Re- 
stricting the time window for the look back but allowing for early exercise 
adds American features to these rather common options. Another virtue of 
the presented Monte Carlo algorithm is to allow option valuation for more 
complicated stochastic processes than Wiener-type diffusion. 
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Table I: Parameters for random ideal 

Wiener processes and put options 



Parameters: 


r 


a 


*So 


X 


T 


mean value 


0.10 


0.40 


100.0 


100.0 


0.50 


standard deviation 


0.05 


0.20 


50.0 


50.0 


0.25 



Caption: Mean values of riskless interest rate r, volatility a, start price So, 
strike price X and time until expiry T are given in first row, their standard 
deviations respectively in second row. The random values are chosen accord- 
ing to independent normal distributions. In a second step, some parameter 
values may be rejected and replaced by other random values if they are non- 
positive. Furthermore, strike price X is required to be within a 2 (1) -ay/T 
range around the mean St value at maturity S + r • T for options contin- 
gent on the current (average) value. This constraint could be relaxed using 
a more sophisticated importance sampling scheme. The number of evenly 
spaced time steps N t for the Monte Carlo simulation and the binomial-tree 
calculations has been fixed arbitrarily to be equal 100. 



26 



Figure Captions: 

Figure 1: 

Monte Carlo estimation of the optimization function O^(S') for a plain 
vanilla put option. Option and stochastic process parameters are chosen 
as in Hull's book (1993), Example 14.1. The Monte Carlo simulation with 
number of time steps equal 50 has been repeated varying the total number 
of sampled paths: 10 2 , 10 3 , 10 4 and 10 5 . Only the function values at time 
step 45 are displayed for 100 values of the argument S'. The location of the 
YES-NO boundary according to a binomial-tree calculation is pointed at by 
an arrow. 

Figure 2: 

Frequency distribution of the relative errors in simulation estimates of a sam- 
ple of plain put option prices. The error is determined from the normalized 
difference to the prices extracted from binomial-tree calculations. The 
Monte Carlo estimates are calculated for the path sample which has been 
used to determine the YES-NO boundary (biased up) as well as for an in- 
dependently generated path sample (biased down) and the average of the 
two samples (MC average). In addition, the figure displays the distribution 
of the difference between upward and downward biased Monte Carlo values 
(lower left panel). Each path sample encompasses 10 5 generated paths per 
option. The Monte Carlo simulation employed the mode in which the 'global 
maximum' of the optimization function Of(S') was searched. 

Figure 3: 

The content of this figure is the same as in the previous figure, with the 
exception of two features. The mode of the Monte Carlo simulation has been 
to look up the maximum of the optimization function on a grid (see text). 
Furthermore, the distribution of errors for the averaged Monte Carlo values 
which looks very similar to the corresponding distributions of upward and 
downward biased estimates is not displayed. Instead, the distribution of the 
errors for the corresponding European option values is displayed (in the lower 
right panel). 
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Figure 4: 

Frequency distribution of deviations between simulation estimates of put 
option prices on geometric (left side) and arithmetic average values (right 
side) and benchmark calculations. The benchmark for the prices on geometric 
averages is provided by the rather accurate binomial-tree method, and for 
the arithmetic averages by the approximate formula eq. (|17|). 
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Figure 1: 
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